#############################################################
# Generalized DiD using Logistic Regression: Crime Concerns
#############################################################

rm(list=ls())

library(Hmisc)
library(readstata13)
library(foreign)
library(lmtest)
library(sandwich)
library(stargazer)
library(msm)
library(ggplot2)
library(did)
set.seed(1111)

#sink("13_generalized_did_concerns_logit.txt")

# functions
vcovCluster <- function(
  model,
  cluster
)
{
  require(sandwich)
  require(lmtest)
  if(nrow(model.matrix(model))!=length(cluster)){
    stop("check your data: cluster variable has different N than model")
  }
  M <- length(unique(cluster))
  N <- length(cluster)           
  K <- model$rank   
  if(M<50){
    warning("Fewer than 50 clusters, variances may be unreliable (could try block bootstrap instead).")
  }
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj  <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
  rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
  return(rcse.cov)
}

round_df <- function(df, digits) {
  nums <- vapply(df, is.numeric, FUN.VALUE = logical(1))
  
  df[,nums] <- round(df[,nums], digits = digits)
  
  (df)
}

############################
# Read and prepare data 
############################

# load individual data
load("final_individual_2023april6.RData")
names(d)

############################
# Concerns binary
############################

# Haiti 
reg1 = glm(d$crime_priority ~ d$t_stagg_haiti + as.factor(d$county) + as.factor(d$survey), family = "binomial")
reg1_c = round(coeftest(reg1, vcov = vcovCluster(reg1, cluster = d$county)),4)
round(reg1_c[1:6,],4)

# Haiti + covariates
reg2 = glm(d$crime_priority ~ d$t_stagg_haiti + d$gender_respondent + d$age_bins_respondent + d$education_respondent + as.factor(d$county) + as.factor(d$survey), family = "binomial")
reg2_c = round(coeftest(reg2, vcov = vcovCluster(reg2, cluster = d$county)),4)
round(reg2_c[1:6,],4)

# All visas 
reg3 = glm(d$crime_priority ~ d$t_stagg_all + + as.factor(d$county) + as.factor(d$survey), family = "binomial")
reg3_c = round(coeftest(reg3, vcov = vcovCluster(reg3, cluster = d$county)),4)
round(reg3_c[1:6,],4)

# All visas + covariates
reg4 = glm(d$crime_priority ~ d$t_stagg_all + d$gender_respondent + d$age_bins_respondent + d$education_respondent + as.factor(d$county) + as.factor(d$survey), family = "binomial")
reg4_c = round(coeftest(reg4, vcov = vcovCluster(reg4, cluster = d$county)),4)
round(reg4_c[1:6,],4)

############################
# Concerns continuous
############################

# Haiti 
reg5 = glm(d$crime_priority ~ d$change_haiti_s + as.factor(d$county) + as.factor(d$survey), family = "binomial")
reg5_c = round(coeftest(reg5, vcov = vcovCluster(reg5, cluster = d$county)),4)
round(reg5_c[1:6,],4)

# Haiti + covariates
reg6 = glm(d$crime_priority ~ d$change_haiti_s + d$gender_respondent + d$age_bins_respondent + d$education_respondent + as.factor(d$county) + as.factor(d$survey), family = "binomial")
reg6_c = round(coeftest(reg6, vcov = vcovCluster(reg6, cluster = d$county)),4)
round(reg6_c[1:6,],4)

# All visas 
reg7 = glm(d$crime_priority ~ d$change_all_s + + as.factor(d$county) + as.factor(d$survey), family = "binomial")
reg7_c = round(coeftest(reg7, vcov = vcovCluster(reg7, cluster = d$county)),4)
round(reg7_c[1:6,],4)

# All visas + covariates
reg8 = glm(d$crime_priority ~ d$change_all_s + d$gender_respondent + d$age_bins_respondent + d$education_respondent + as.factor(d$county) + as.factor(d$survey), family = "binomial")
reg8_c = round(coeftest(reg8, vcov = vcovCluster(reg8, cluster = d$county)),4)
round(reg8_c[1:6,],4)

###############################
# Tables
###############################

# Table a5
stargazer(reg5_c,reg6_c,reg7_c,reg8_c,
          digits = 3,
          title="",
          align=TRUE, 
          omit.stat=c("LL","ser","f","rsq","adj.rsq"), 
          no.space=TRUE,  
          multicolumn = TRUE,
          table.placement = "H",
          column.separate = c(2,2),
          covariate.labels=c("Visas from Haiti","Visas from all countries"),
          dep.var.caption  = "Crime Concerns" ,
          dep.var.labels.include = FALSE,
          star.cutoffs = c(0.05,0.01,0.001),
          omit = c("Constant","gender_respondent","age_bins_respondent","education_respondent","county","survey"),
          add.lines = list(c("Covariates","No","Yes","No","Yes"),c("Observations","13,551","13,551","13,551","13,551")),
          out="tableA5.txt")

# Table a5
stargazer(reg1_c,reg2_c,reg3_c,reg4_c,
          digits = 3,
          title="",
          align=TRUE, 
          omit.stat=c("LL","ser","f","rsq","adj.rsq"), 
          no.space=TRUE,  
          multicolumn = TRUE,
          table.placement = "H",
          column.separate = c(2,2),
          covariate.labels=c("Visas from Haiti","Visas from all countries"),
          dep.var.caption  = "Crime Concerns" ,
          dep.var.labels.include = FALSE,
          star.cutoffs = c(0.05,0.01,0.001),
          omit = c("Constant","gender_respondent","age_bins_respondent","education_respondent","county","survey"),
          add.lines = list(c("Covariates","No","Yes","No","Yes"),c("Observations","13,551","13,551","13,551","13,551")),
          out="tableA6.txt")

sink()

